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Residence Time Distribution of Sand Grains in the 
1 -Dimensional Abelian Sandpile Model 



Punyabrata Pradhan and Apoorva Nagar 



Abstract — We study the probability distribution of residence time, T, of 
the sand grains in the one dimensional abelian sandpile model on a lattice 
of L sites, for T « L 2 and T » L 2 . The distribution function decays 
as exp(— )• We numerically calculate the coefficient Kt, for the value 
of L upto 150 . Interestingly the distribution function has a scaling form 
T^rfi Jfi ) w '*h a ^ b f° r large L. 
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I. Introduction 

There are many slowly driven extended dissipative systems 
in nature which organize themselves to a critical state. These 
are called self organized critical (SOC) systems. Such systems 
show long-ranged spatial and temporal correlations, i.e., power 
law behaviour in the correlation functions; for example, the fre- 
quency of earthquakes as a function of energy released, the mean 
square velocity difference as a function of distance in the tur- 
bulant fluids, height-height correlation in the case of surface 
roughening of growing interfaces 0, 0, 0, 0- The time 
series of electrical noise shows power law in the power spec- 
trum ( "1 //" noise) [2 1. Sandpile models are the simplest mod- 
els where important features of Self Organized Critical (SOC) 
states are present. Abelian sandpile model is a subclass of these. 
There are many papers |3|, |5|, [6|, |7| concerned mainly in 
studying the distributions of avalanche sizes and avalanche du- 
rations and the corresponding critical exponents, but there is not 
much work in the literature where the distribution of times these 
grains spend inside the system has been studied theoretically. 
Although it has been studied in an experiment on ricepiles and in 
simulation in one dimension [6 1. In this paper we study the resi- 
dence time distribution of a marked sand grain, when the system 
is driven by adding sand grains at one end (say, left). Our model 
is BTW sandpile model in 1 -dimension which is simplest 
to study. We find out the scaling form to be -p-J(-p-) for large 
L, where T is the residence time and L is the size of the sys- 
tem. This scaling function f(x), although it diverges as x — > 0, 
is consistent with the normalization condition, since lower limit 
of x has a cutoff at a; = -p-. The residence time distribution 
decays exponentially for large value of T. We determine nu- 
merically the coefficient, Kl, of TJ L 2 in the exponential decay 
for some finite sizes of the 1-diemnsional lattice. 



II. Model 

We take an one dimensional lattice of L sites where sites are 
indexed by i with 1 < i < L. The configuration of the sand- 
pile is specified by height variable (number of sand grains) Zj at 
any site i. The dynamics is as follows. Whenever the height is 
higher than some critical value z c , the site topples and two sand 



grains go to the nearest neighbour sites, one goes to the left and 
the other goes to the right. Which one goes to the left is chosen 
at random. Here we consider both z c = 1 and z c — 2. There 
is one layer of sand grains for z c = 1. For z c = 2, there are 
two layers of grains. When the height of the end site at i = 1 or 
i = L crosses the critical height, the site topples. One sand grain 
goes out of the system and other goes to the next neighbour site. 
No sand is added till the system comes to the stable configura- 
tion, i.e., there is no more topplings in the system. The time unit 
is measured as the duration between the successive addition of 
two sand grains. The residence time of a sand grain is defined 
as the time spent by that grain inside the system. P(T\L) is the 
probability distribution function of the residence time T where 
L is the number of sites in the 1-dimensoinal lattice. 
First we consider the model with critical height z c = 1 and start 
from the configuration with the height z = at all sites. Now 
we keep on adding sand grain at the leftmost site. After first ad- 
dition of grain, there is no toppling. Then we add the next sand 
grain. Now there will be toppling, one sand grain will leave the 
system, other grain will move to the right and stays there until 
next grain comes to the site, and so on. So the evolution rule is 
totally deterministic. After some time (~ L 2 ) the system falls 
into a cycle. This cycle, consisting of L + 1 stable configu- 
rations, has at most one site with = and other sites with 
Zi — 1 (recurrent configuration space). In the steady state of the 
system, all the configurations are equally probable. 
Of the L + 1 recurrent configurations, L configurations have 
a site with zero height and one configuration has all sites with 
height 1. If there is a configuration which has i-th site with zero 
height, in the next time step the configuration changes to one 
which has (i — l)-th site with height zero (If i = 1, in the next 
step all heights become 1). At the next time step, the (i — 2)-th 
site is with height zero. After some time there is 1st site with 
zero height. When next sand grain is added, all sites become 
with height 1. In the next step, last site becomes with height 
zero and so on. Even though time evolution of the system is 
completely deterministic, movement of sand grains is stochas- 
tic. Whenever any site topples, one of the two sand grains goes 
to the left nearest neighbour and stops, and other grain moves to 
the right nearest neighbour. The probability that one particular 
grain of the two goes to the left and other goes to the right is 
1/2. If the height at the right site is zero, toppling stops, oth- 
erwise, there is another toppling. If any sand grain falls on the 
site with height zero, it stops moving and waits till the site with 
height zero is at the left of it. When the L-th site has the height 
zero, the sand grain starts moving again in the similar fashion. 
The waiting times are not same at every site and waiting times 
are correlated. 
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III. Analytical Results 

First we shall discuss the case for z c = 1. When the resi- 
dence time T is very small compared to L, sand grains always 
get ejected from the left and do not go beyond the site i = T 
(except for the case when Zi = 1 at all sites at time T and grain 
might be ejected from the right end). It's very unlikely that the 
marked sand grain meets the site with height zero (the probabil- 
ity — (L+i) wri i c h S oes to zero f° r large value of L). In the limit 
of L large compared to T, the residence time distribution is well 
approximated by the distribution of the first passage (at i = 1) 
time of a simple unbiased random walker with 2T steps random 
walk. So, for T << L, We can exactly calculate the residence 
time distribution 1 1 1 1 



P(T\L) = 



2T\ 



T!(T+ 1)! 



2 -(2T+l) 



where P(T\L) is the probability that the marked sand grain gets 
ejected from the system immediatly after time T. Using Ster- 
ling approximation to the above expression, we get the function 
P(T\L) equals to ^=5^72 for 1 << T << L. Now we can 
extend our result whenever T << L 2 . The argument is as fol- 
lows. After time steps T ~ L a where a < 2, the standard 
deviation of the position of the marked sand grain goes as L a l 2 . 
So the probability that the marked grain will meet the site with 
height zero goes as L Q / 2_1 which will tend to zero in the large 
L limit. Therefore for any 1 << T « L 2 , the limiting proba- 
bility distribution P(T\L) is proportional to T~ 3 / 2 , the first pas- 
sage (at i = l) time distribution of a simple unbiased random 
walker with 2T steps random walk. For the large value of T the 
distribution function decays exponentially like exp(—K lT/L 2 ) 
where Kl depends on lattice size. Later we will describe how 
to compute the coefficient Kl for any finite L. 

Even though we don't know the distribution function for all 
T, we can find out the first moment of the residence time distri- 
bution easily. We define mass of the sandpile as the total number 
of particles in the pile (i.e., X)i=i ^'s easy to see that mean 
residence time (T) = (Total mass of the pile). To prove this, let 
us define an indicator funtion rj n ,T as given below. 

Vn,t — 1 if the sand grain, added at time n, is in the system 

at time t, otherwise rj n t = 0. 
The mean residence time can be written as 

1 N T T M 

{T) = aa.^oo n £ £ = m £ £ ^ 

71=1 t=l t=l 71=1 

! r 

= -jj 2_j (Total mass of the pile at time t) 

= ( Total mass of the pile ) 

The average of the total mass in the pile is [j^p[(L — 1) + j^p[L] 
which goes as L for large L. 

We find out the scaling form of the function P(T\L) using the 



mean residence time and the previous limiting distribution. The 
scaling finction f(x) is defined as 

f(x)dx = lim i a " 6 Prob(.T J L b <T <(x + dx)L b ) 

L — >oo 

So the form of the distribution function P(T\L) is -^sfijp)- 
Since P(T\L) varies as T~ 3 / 2 for 1 << T << L, f(x) must 
goes as ir~ 3 / 2 for very small value of x and therefore 77- /(-p-) 
goes as £(3b/2-a) r -3/2. As P(T\L) is independent of L for 
1 << T « L, we get a/b = 3/2. Even though the scaling 
function f(x) is divergent as x tends to zero, we can normalize 
P(T\L) without any problem, since x has lower cutoff -ig. As 
f(x) goes as a; -3 / 2 , the normalisation integral is 



L 



1/L b 



f{x)dx - i 36 / 2 



which is independent of L as a/b = 3/2. The mean residemce 
time (T) is given by 



f(T/L b )TdT - I}' 2 



L 3b/2 



As there is no divergence in the integrand, we can put the 
lower limit in the integration zero. Since (T) ~ L, we find 
6 = 2 and therefore can write the scaling form of the P(T\L) 
as T~ 3 / 2 /(-p-) for large L. In Fig.l. probability distributions, 
P(T\L), of residence time of sand grains are plotted for differ- 
ent lattice sizes L = 25,50, 100, 150. In Fig. 2. all the curves 
for P(T\L) collapse to a single curve when we plot L 3 P(T\L) 
against T/L 2 . 

IV. Computational Results 

Now we study the behaviour of the function P(T\L) when 
T >> L 2 . Since time evolution of states of the system is a 
Markov process, there is exponential decay of probabilities of 
states for large value of time, due to the presence of the absorb- 
ing boundaries. So the probability distribution P(T\L) must 
decays as exp(— K^T/ L 2 ) for very large T. For large L, the co- 
efficient, Kl, tends to a constant K and scaling function f(x) 
must goes as exp(— Kx) for large x. We define K as a limit 
given below, 



K = - lim 

L — >oo 



lim 



L 2 \nP{T\L) 
f 



We numerically calculate values of Kl 's for various finite lattice 
sizes by defining transition probability matrices for the system 
going from one recurrent configuration to another. Now we need 
to distinguish all configurations with a specific site which has 
height zero. Configurations are distinguishable with respect to 
the position of the marked grain. Marked grain can be at any site 
except the site where height is zero. For simplicity, we construct 
the transition matrices for small value of L, say L = 4. When 
the site with zero height is at the end, i.e., at i = 4, three distin- 
guishable configurations are represented respectively as | * 1 10), 
|1 * 10), |11 * 0) where "Star" denotes the position of marked 
grain, "1" denotes the site with height one and "0" denotes the 
site with height zero. If we keec on adding sand at the left. 
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the transitions will occur from one state to another. In this case 
transitions will be from 

{1110} -> {1101} -> {1011} -» {0111} 

-» {1111} -» {mo} 

(same initial states after L + 1 = 5 time steps). 

Here marked grain can be at any one of the sites with height 1. 
As mentioned before, we represent basis states for a particular 
configuration as \j, Y ) where j denotes the position of marked 
particle and Y denotes the index of the site with height zero. 
Since marked grain cannot stay at the site with height zero, there 
are L — 1 basis states for a particular configuration (i.e. for 
a fixed Y ). Y can take value from to L. Y = means 
that all sites are with height 1. Whenever we add sand grain 
at the left end, there is transition from any of the basis states 
with some value of Y to any of the basis states with Yq — 1 (i.e., 
\j,Y ) — * \i, Yq — 1)). i and j denote the position of marked sand 
grain. Now we define elements of the transitoin matrix,^,, (i\ j), 
as the transition probability from the j-th state to z-th state where 
Y denotes that the transition occurs from the configuration with 
height zero at Yo-th site to the configuration with height zero at 
(Y — l)-th site. The transition matrices %, %, T 2 , T\, % for 
L = 4 are written below explicitly. 
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The transition matrix = YIy =o ^~ y q gi yes back the ini- 
tial configuration after the period of time L + 1. The super- 
script in 7~( L ) is denoting that the marked grain is added when 
the initial configuration is with Yq = L. In general, if the 
grain is added with the initial configuration where position of 
the site with heigt zero is at z = Yq, the transition matrix will 
beTW = T Yo _ 1 T Yq _ 2 ....T l _ 1 T l T .. ..T Yo -iT Yq . The dif- 
ferent values of Y in transition matrices correspod to the addi- 
tion of sand grains to the configurations with different sites with 
height zero. To find out the coefficient Kl in the exponential, 
we diagonalize any matrix T^ Yo \ say T^ L \ and take the largest 
eigenvalue, X max - We diagonalise any one matrix of the differ- 
ent transitoin matrices because the eigen values for the differ- 
ent transitoin matrices, T^ Yo h (we get various T^ Y °h just with 
cyclic permutation of matrices %, T 2 , %....Tl etc.) are same. 
For large value of time T (>> L 2 ), the residence time distribu- 
tion decays exponentially like (X max ) T ^ L+1 ^ which equals to 
exp( I ^ T In X max )- The coefficient K L is then -j^j In X max . 
In the Fi2.3. we have plotted Kr aeainst the number of sites in 




Residence Time T 
Fig. 1 

THE PROBABILITY DISTRIBUTION OF RESIDENCE TIME 
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Fig. 2 

The plot of scaling function f(x) showing collapse of all the 

CURVES FOR DIFFERNET LATTICE SIZES TO A SINGLE CURVE. 



the 1 -dimensional lattice, L. It shows that K L saturates to the 
value of K. So the coefficient is independent of lattice size L 
for large value of L. 



We can easily extend all our previous result for z c = 1 to 
z c = 2. The toppling rules for z c = 2 are similar to the pre- 
vious case z c = 1 except that sand grain can stay at the site 
even if there is a toppling. Whenever height at any site is three, 
the site gets unstable and topples. One sand grain of the three 
goes to right nearest neighbour site, one goes to the left nearest 
neighbour site and the other one stays put. It is decided at ran- 
dom which one goes to the right, which one goes to the left. The 
probablity that a particular grain of the three going to the right or 
to the left, or staying there is 1 /3. The height of any site is either 
1 or 2, when the system is stable. In the steady state, at most one 
site has height 1 and all the others have height 2. All configura- 
tions in the recurrent configuration soace are eauallv crobable. 
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Fig. 3 

Plot of K l versus lattice size L 



The number of recurrent configuratins is same as before. The 
residence time distribution P(T\L), for 1 << T « L, is 

T 

°2TO T!(T+1)! T 3 / 2 

(1) 

Here also we can extend this result for 1 << T << L 2 , just 
as we did in the case for z c = 1. The probability distribution 
functin has the same scaling form as for z c = 1. 

To calculate the coefficient Kl, we can construct the trasition 
probability matrix exactly in the similar way for the previous 
case. In this case the number of basis states is L and Kl can be 
determined as before. 



V. Summary 

In summary, we have calculated probability distribution of the 
residence time T in two limitting cases T « L 2 and T > > L 2 , 
have shown that the distribution function has a scaling form 
-p-/(-p-) and computed the coefficient of T/L 2 , Kl, in the ex- 
ponential decay of the distribution function for the values of L 
upto 150. 
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